Ce document présente un workflow complet pour la délimitation de bassins versants à partir d’un Modèle Numérique d’Élévation (MNE). Le processus inclut le prétraitement du MNE, la génération de réseaux hydrographiques et la délimitation finale des bassins versants.
Chargement des packages
#Packages de base pour la manipulation de données spatialesoptions(repos =c(CRAN ="https://cloud.r-project.org"))install.packages("Rcpp")
The downloaded binary packages are in
/var/folders/qb/stf0vsnn0zn0509g0d3ch4r40000gn/T//Rtmp04WMM7/downloaded_packages
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(raster)
Loading required package: sp
Attaching package: 'raster'
The following object is masked from 'package:dplyr':
select
library(terra) # Alternative moderne au package raster
terra 1.8.21
Attaching package: 'terra'
The following object is masked from 'package:tidyr':
extract
library(sf) # Pour les objets vectoriels
Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(tmap) # Cartographie interactivelibrary(stars) # Autre format de données spatiales
Loading required package: abind
library(whitebox) # Interface avec WhiteboxTools# Initialisation de WhiteboxToolswhitebox::wbt_init()wbt_version()
[1] "WhiteboxTools v1.4.0 by Dr. John B. Lindsay (c) 2017-2020"
[2] ""
[3] "WhiteboxTools is an advanced geospatial data analysis platform developed at"
[4] "the University of Guelph's Geomorphometry and Hydrogeomatics Research "
[5] "Group (GHRG). See https://jblindsay.github.io/ghrg/WhiteboxTools/index.html"
[6] "for more details."
# Configuration des thèmes et modes d'affichagetheme_set(theme_classic())tmap_mode("view")
ℹ tmap mode set to "view".
Importation et préparation du MNE
Cette section importe le Modèle Numérique d’Élévation (MNE) et effectue les premiers traitements pour éliminer les valeurs aberrantes.
# Lecture du fichier MNE (format .adf)MNE <-rast("/Users/thejo/Desktop/TD/Donnees_Audrey/DEM/demwatershed")# Nettoyage des valeurs négatives (souvent des artefacts)MNE[MNE <0] <-NA
[v3->v4] `tm_raster()`: instead of `style = "cont"`, use col.scale =
`tm_scale_continuous()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values') to
'tm_scale_continuous(<HERE>)'
! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.
SpatRaster object downsampled to 3740 by 2675 cells.
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "PuOr" is named
"brewer.pu_or"
2. Conversion et ombrage du MNE
Conversion au format .tif pour une meilleure compatibilité et génération d’un ombrage pour la visualisation.
# Conversion en format .tifwriteRaster(MNE, "/Users/thejo/Desktop/TD/trois/donnees/MNE.tif", overwrite =TRUE)
[v3->v4] `tm_raster()`: instead of `style = "cont"`, use col.scale =
`tm_scale_continuous()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values') to
'tm_scale_continuous(<HERE>)'
[v3->v4] `tm_raster()`: use `col.legend = tm_legend_hide()` instead of
`legend.show = FALSE`.
! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.
SpatRaster object downsampled to 3740 by 2675 cells.
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "-Greys" is named
"greys" (in long format "brewer.greys")
Correction des dépressions
Les dépressions dans le MNE peuvent affecter l’analyse hydrologique. Cette section les corrige.
# Conversion en format float32 avec valeur NA spécifiqueterra::writeRaster( MNE, "/Users/thejo/Desktop/TD/trois/donnees/MNE_float32.tif", datatype ="FLT4S", overwrite =TRUE, NAflag =-9999)
# Correction des dépressions par la méthode du moindre coûtwbt_breach_depressions_least_cost(dem ="/Users/thejo/Desktop/TD/trois/donnees/MNE_float32.tif",output ="/Users/thejo/Desktop/TD/trois/donnees/MNE_breach_dep.tif",dist =100, # Distance maximale de correction, à ajuster selon le casfill =TRUE)
[1] "breach_depressions_least_cost - Elapsed Time (excluding I/O): 34min 54.999s"
# Remplissage complémentaire des dépressionswbt_fill_depressions_wang_and_liu(dem ="/Users/thejo/Desktop/TD/trois/donnees/MNE_breach_dep.tif",output ="/Users/thejo/Desktop/TD/trois/donnees/MNE_cleaned.tif")
[1] "fill_depressions_wang_and_liu - Elapsed Time (excluding I/O): 1min 22.107s"
MNE_clean <-rast("/Users/thejo/Desktop/TD/trois/donnees/MNE_cleaned.tif")# Vérification des résultatssummary(values(MNE_clean))
MNE_cleaned
Min. : 0
1st Qu.:475
Median :541
Mean :500
3rd Qu.:591
Max. :944
NA's :131419502
Analyse du réseau hydrographique
Calcul de l’accumulation de flux et de la direction d’écoulement.
# Accumulation de flux (méthode D8)wbt_d8_flow_accumulation(input ="/Users/thejo/Desktop/TD/trois/donnees/MNE_cleaned.tif",output ="/Users/thejo/Desktop/TD/trois/donnees/flow_accD8.tif")
[1] "d8_flow_accumulation - Elapsed Time (excluding I/O): 11.71s"
# Transformation logarithmique pour meilleure visualisation sur ArcGIS proflow_acc <-rast("/Users/thejo/Desktop/TD/trois/donnees/flow_accD8.tif")scaled_acc <-log1p(flow_acc)
# Direction d'écoulement (méthode D8)wbt_d8_pointer(dem ="/Users/thejo/Desktop/TD/trois/donnees/MNE_cleaned.tif",output ="/Users/thejo/Desktop/TD/trois/donnees/pointer3.tif")
[1] "d8_pointer - Elapsed Time (excluding I/O): 3.470s"
pointer3 <-rast("/Users/thejo/Desktop/TD/trois/donnees/pointer3.tif")# Vérifier les valeurs du flow direction (Ne devrait que contenir: 1,2,4,8,16,32,64,128)freq(pointer3)
Création et préparation des points d’exutoire pour la délimitation des bassins versants.
# Création des points d'exutoireppoints <-tribble(~Lon, ~Lat,-63.369243, 50.828091)# Conversion en objet spatialppointsSP <-SpatialPoints(ppoints, proj4string =CRS("+proj=longlat +datum=WGS84"))# Export en shapefileshapefile(ppointsSP, filename ="pourpoints.shp", overwrite =TRUE)
Extraction du réseau hydrographique
Seuillage de l’accumulation de flux pour identifier les cours d’eau.
wbt_extract_streams(flow_accum ="/Users/thejo/Desktop/TD/trois/donnees/flow_accD8.tif",output ="/Users/thejo/Desktop/TD/trois/donnees/raster_streams.tif",threshold =200# Ajuster selon la résolution et la zone d'étude)
[1] "extract_streams - Elapsed Time (excluding I/O): 2.392s"
wbt_jenson_snap_pour_points(pour_pts ="/Users/thejo/Desktop/TD/pourpoints.shp",streams ="/Users/thejo/Desktop/TD/trois/donnees/raster_streams.tif",output ="/Users/thejo/Desktop/TD/trois/donnees/snappedpp.shp",snap_dist =200# Distance d'ajustement (unités du SCR))
[1] "jenson_snap_pour_points - Elapsed Time (excluding I/O): 0.2s"
snappp <-shapefile("/Users/thejo/Desktop/TD/trois/donnees/snappedpp.shp")streams <-raster("/Users/thejo/Desktop/TD/trois/donnees/raster_streams.tif")#Vérifier concordance des systèmes de coordonnéesprint("Pointer3 CRS:")
[1] "Pointer3 CRS:"
print(crs(pointer3))
[1] "PROJCRS[\"unnamed\",\n BASEGEOGCRS[\"Unknown datum based upon the GRS 1980 ellipsoid\",\n DATUM[\"Not specified (based on GRS 1980 ellipsoid)\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101004,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4019]],\n CONVERSION[\"Lambert Conic Conformal (2SP)\",\n METHOD[\"Lambert Conic Conformal (2SP)\",\n ID[\"EPSG\",9802]],\n PARAMETER[\"Latitude of false origin\",49,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8821]],\n PARAMETER[\"Longitude of false origin\",-95,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8822]],\n PARAMETER[\"Latitude of 1st standard parallel\",49,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8823]],\n PARAMETER[\"Latitude of 2nd standard parallel\",77,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8824]],\n PARAMETER[\"Easting at false origin\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8826]],\n PARAMETER[\"Northing at false origin\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8827]]],\n CS[Cartesian,2],\n AXIS[\"easting\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"northing\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]]"
[v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
visual variable `col` namely 'palette' (rename to 'values') to col.scale =
tm_scale(<HERE>).
SpatRaster object downsampled to 3740 by 2675 cells.
The visual variable "col" of the layer "raster" contains a unique value. Therefore a discrete scale is applied (tm_scale_discrete).
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
"brewer.blues"